#SIMULATIONS FOR: ENVIRONMENTAL POLICY IN GENERAL EQUILIBRIUM: NEW INSIGHTS FROM A CANONICAL MODEL

#PIERRE MEREL, COMPLETED JUNE 2021 FOR JAERE

#This code seeks to replicate results in Table 2 and Table 3 of Fullerton & Heutel (F&H 2007 JPubE) and generate the tables 
#in Garnache and Merel (JAERE). We could not exactly replicate the results of F&H,
#likely due to a coding error in their equilibrium system denominator, but we get close. Our numerical results 
#match their analytical formulas (and ours). We use our numerical results for comparisons.
#The code also provides support for various claims made throughout our paper and its online appendix.

library(xtable)


#Allen elasticity functions for own-price effects
ell<-function(fk,fz,tl,tk,tz){
  e<--1/tl*(tk*fk+tz*fz)
  return(e)
}
ekk<-function(fl,fz,tl,tk,tz){
  e<--1/tk*(tl*fl+tz*fz)
  return(e)
}
ezz<-function(fl,fk,tl,tk,tz){
  e<--1/tz*(tl*fl+tk*fk)
  return(e)
}


#######################################MAIN ANALYSIS###################################

##TABLE 2 IN F&H 2007 (SECTION 5)

eLK<-1
sX<-1
sU<-1
gL<-0.25
gK<-0.25
tZ<-0.25
tL<-0.45
tK<-1-tL-tZ
sL<-gK*tL/(gK*tL+gL*tK)
sK<-1-sL
d<-gL*sL/tL
Y<-d/(1+d)
L<-sL*(1+gL)/(sL*(1+gL)+sK*(1+gK))
TZHAT<-0.10

#The function Table2 generates all tables that investigate alternative choices of Allen substitution elasticities
#also computes eKK, eZZ, and subdeterminant to rule out choices that violate negative semidefiniteness of the Slutsy matrix
#also computes C1 C2 C3 to determine drivers of discrepancies across numeraire choices


Table2<-function(ekz,elz,a,b) {
  eLL<-ell(eLK,elz,tL,tK,tZ)
  eKK<-ekk(eLK,ekz,tL,tK,tZ)
  eZZ<-ezz(elz,ekz,tL,tK,tZ)
  delta<-sX*(1-b*tZ)*(1+gL*sL+gK*sK)+sU*(gL-gK)*(tL-sL+tZ*(b*sL+(1-b)*a))-eLL*tL*gL*(1+gK)*(b*tK+(1-b)*(1-a))-eKK*tK*gK*(1+gL)*(b*tL+(1-b)*a)+eLK*(tL*gK*(1+gL)*(b*tK+(1-b)*(1-a))+tK*gL*(1+gK)*(b*tL+(1-b)*a))
  d2<-tZ*tK/sL*(eKK*gK+eZZ*(gK-d)^2/gK-2*ekz*(gK-d))
  determinant<-eKK*eZZ-(ekz)^2
  C1<-(1-b*tZ)*(1+gL*sL+gK*sK)
  C2<-(gL-gK)*(tL-sL+tZ*(b*sL+(1-b)*a))
  C3<--eLL*tL*gL*(1+gK)*(b*tK+(1-b)*(1-a))-eKK*tK*gK*(1+gL)*(b*tL+(1-b)*a)+eLK*(tL*gK*(1+gL)*(b*tK+(1-b)*(1-a))+tK*gL*(1+gK)*(b*tL+(1-b)*a))
  pollu<-100*TZHAT/delta*((1+d)*tZ*eZZ*sX+d2*sU-tZ*sX*sU-(gL*tK+gK*tL+gL*gK)*(tL*eLK*elz+tK*eLK*ekz+tZ*elz*ekz))
  source<-100*tZ/delta*(sU*(gK-gL)+gL*(1+gK)*elz-gK*(1+gL)*ekz)*TZHAT
  use<-100*tZ/delta*(sX*(1+gL*sL+gK*sK)+(tL*gK*(1+gL)+tK*gL*(1+gK))*eLK+gL*(1+gK)*(sK-tK)*elz+gK*(1+gL)*(sL-tL)*ekz)*TZHAT
  return(cbind(pollu,source,use,C1,C2,C3,determinant))
}

M2<-data.frame(ekz=c(1,0,.5,1,-0.5,0,0.5,1,-0.5,0,0.5,1),
                 elz=c(-0.5,0,0,0,0.5,0.5,0.5,0.5,1,1,1,1))
M2$pX<-Table2(M2$ekz,M2$elz,sL,0)
M2$pY<-Table2(M2$ekz,M2$elz,1,1)
M2$merel<-Table2(M2$ekz,M2$elz,sL,Y)
M2$w<-Table2(M2$ekz,M2$elz,1,0)
M2$r<-Table2(M2$ekz,M2$elz,0,0)
M2$merel2<-Table2(M2$ekz,M2$elz,L,0)
M2$mean<-1/6*(M2$pX+M2$pY+M2$merel+M2$w+M2$r+M2$merel2)
M2$sd<-sqrt(1/6*((M2$pX-M2$mean)^2+(M2$pY-M2$mean)^2+(M2$merel-M2$mean)^2+(M2$w-M2$mean)^2+(M2$r-M2$mean)^2+(M2$merel2-M2$mean)^2))

# Convert to latex format
xtable(cbind(M2$ekz,M2$elz,M2$pX[,1],M2$pY[,1],M2$merel[,1],M2$w[,1],M2$r[,1],M2$sd[,1]),caption = "pollu",digits=c(1,1,1,2,2,2,2,2,2))
xtable(cbind(M2$ekz,M2$elz,M2$pX[,2],M2$pY[,2],M2$merel[,2],M2$w[,2],M2$r[,2],M2$sd[,2],M2$pX[,3],M2$pY[,3],M2$merel[,3],M2$w[,3],M2$r[,3],M2$sd[,3]),caption = "sourceANDuse",digits=c(1,1,1,2,2,2,2,2,2,2,2,2,2,2,2))
#NOTE: Remember to eliminate rows for which the subdeterminant is negative.


##TABLE THAT COMPUTES THE TAX RATE EQUIVALENT TO A 10% INCREASE IN THE TAX RELATIVE TO THE CPI (APPENDIX D.1)

TableTaxRate<-function(ekz,elz,a,b) {
  eLL<-ell(eLK,elz,tL,tK,tZ)
  eKK<-ekk(eLK,ekz,tL,tK,tZ)
  eZZ<-ezz(elz,ekz,tL,tK,tZ)
  aa<-sL
  bb<-Y
  deltabenchmark<-sX*(1-bb*tZ)*(1+gL*sL+gK*sK)+sU*(gL-gK)*(tL-sL+tZ*(bb*sL+(1-bb)*aa))-eLL*tL*gL*(1+gK)*(bb*tK+(1-bb)*(1-aa))-eKK*tK*gK*(1+gL)*(bb*tL+(1-bb)*aa)+eLK*(tL*gK*(1+gL)*(bb*tK+(1-bb)*(1-aa))+tK*gL*(1+gK)*(bb*tL+(1-bb)*aa))
  delta<-sX*(1-b*tZ)*(1+gL*sL+gK*sK)+sU*(gL-gK)*(tL-sL+tZ*(b*sL+(1-b)*a))-eLL*tL*gL*(1+gK)*(b*tK+(1-b)*(1-a))-eKK*tK*gK*(1+gL)*(b*tL+(1-b)*a)+eLK*(tL*gK*(1+gL)*(b*tK+(1-b)*(1-a))+tK*gL*(1+gK)*(b*tL+(1-b)*a))
  THAT<-100*TZHAT*delta/deltabenchmark
  return(THAT)
}

MTAX<-data.frame(ekz=c(1,0,.5,1,-0.5,0,0.5,1,-0.5,0,0.5,1),
               elz=c(-0.5,0,0,0,0.5,0.5,0.5,0.5,1,1,1,1))
MTAX$pX<-TableTaxRate(MTAX$ekz,MTAX$elz,sL,0)
MTAX$pY<-TableTaxRate(MTAX$ekz,MTAX$elz,1,1)
MTAX$merel<-TableTaxRate(MTAX$ekz,MTAX$elz,sL,Y)
MTAX$w<-TableTaxRate(MTAX$ekz,MTAX$elz,1,0)
MTAX$r<-TableTaxRate(MTAX$ekz,MTAX$elz,0,0)
MTAX$merel2<-TableTaxRate(MTAX$ekz,MTAX$elz,L,0)

xtable(cbind(MTAX$ekz,MTAX$elz,MTAX$pX,MTAX$pY,MTAX$merel,MTAX$w,MTAX$r),caption = "equivalent tax",digits=c(1,1,1,2,2,2,2,2))


##TABLE 3 IN F&H 2007  (APPENDIX D.2)

#The function Table3 generates a table that investigates alternative choices of gammaK-gammaL
#also computes eKK, eZZ, and subdeterminant to rule out choices that violate NSD of the Slutsky matrix

rm(gL,gK,eLK,tL,tK,tZ,sL,sK,sX,sU,d,Y,L,TZHAT)

TZHAT<-0.10

Table3<-function(g,model){
  ky<-ifelse(!g==0,(-1+sqrt(1+8*g*(2*g+1)/25))/(2*g),0.08)
  kx<-2/5-ky
  ly<-1/5-ky
  lx<-2/5+ky
  gL<-ly/lx
  gK<-ky/kx
  eLZ<-1
  eKZ<--0.5
  eLK<-1
  sX<-1
  sU<-1
  tZ<-0.25
  tL<-ly*15/4
  tK<-ky*15/4
  sL<-lx*5/4
  sK<-1-sL
  d<-gL*sL/tL
  Y<-d/(1+d)
  L<-sL*(1+gL)/(sL*(1+gL)+sK*(1+gK))
  
  eLL<-ell(eLK,eLZ,tL,tK,tZ)
  eKK<-ekk(eLK,eKZ,tL,tK,tZ)
  eZZ<-ezz(eLZ,eKZ,tL,tK,tZ)
  if (model=="pX") {b<-0
  a<-sL}
  if (model=="pY") {b<-1
  a<-1}
  if (model=="w") {b<-0
  a<-1}
  if (model=="r") {b<-0
  a<-0}
  if (model=="merel") {b<-Y
  a<-sL}
  if (model=="merel2") {b<-0
  a<-L}
  d2<-tZ*tK/sL*(eKK*gK+eZZ*(gK-d)^2/gK-2*eKZ*(gK-d))
  delta<-sX*(1-b*tZ)*(1+gL*sL+gK*sK)+sU*(gL-gK)*(tL-sL+tZ*(b*sL+(1-b)*a))-eLL*tL*gL*(1+gK)*(b*tK+(1-b)*(1-a))-eKK*tK*gK*(1+gL)*(b*tL+(1-b)*a)+eLK*(tL*gK*(1+gL)*(b*tK+(1-b)*(1-a))+tK*gL*(1+gK)*(b*tL+(1-b)*a))
  determinant<-eKK*eZZ-(eKZ)^2
  source<-100*tZ/delta*(sU*(gK-gL)+gL*(1+gK)*eLZ-gK*(1+gL)*eKZ)*TZHAT
  use<-100*tZ/delta*(sX*(1+gL*sL+gK*sK)+(tL*gK*(1+gL)+tK*gL*(1+gK))*eLK+gL*(1+gK)*(sK-tK)*eLZ+gK*(1+gL)*(sL-tL)*eKZ)*TZHAT
  pollu<-100*TZHAT/delta*((1+d)*tZ*eZZ*sX+d2*sU-tZ*sX*sU-(gL*tK+gK*tL+gL*gK)*(tL*eLK*eLZ+tK*eLK*eKZ+tZ*eLZ*eKZ))
  return(cbind(pollu,source,use,tK,sK,eKK,eZZ,determinant))
}

M3<-data.frame(gkgl=c(-0.25,-0.20,-0.15,-0.1,-0.05,0,0.05,0.1,0.15,0.2,0.25))
M3$pX<-Table3(M3$gkgl,"pX")
M3$pY<-Table3(M3$gkgl,"pY")
M3$w<-Table3(M3$gkgl,"w")
M3$r<-Table3(M3$gkgl,"r")
M3$merel<-Table3(M3$gkgl,"merel")
M3$merel2<-Table3(M3$gkgl,"merel2")
M3$mean<-1/6*(M3$pX+M3$pY+M3$merel+M3$w+M3$r+M3$merel2)
M3$sd<-sqrt(1/6*((M3$pX-M3$mean)^2+(M3$pY-M3$mean)^2+(M3$merel-M3$mean)^2+(M3$w-M3$mean)^2+(M3$r-M3$mean)^2+(M3$merel2-M3$mean)^2))

# Convert to latex format
xtable(cbind(M3$gkgl,M3$pX[,4],M3$pX[,5],M3$pX[,1],M3$pY[,1],M3$merel[,1],M3$w[,1],M3$r[,1],M3$merel2[,1],M3$sd[,1]),caption = "M3pollu",digits=c(1,2,4,4,2,2,2,2,2,2,2))
xtable(cbind(M3$gkgl,M3$pX[,2],M3$pY[,2],M3$merel[,2],M3$w[,2],M3$r[,2],M3$merel2[,2],M3$sd[,2],M3$pX[,3],M3$pY[,3],M3$merel[,3],M3$w[,3],M3$r[,3],M3$merel2[,3],M3$sd[,3]),caption = "sourceANDuse",digits=c(1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2))
#NOTE: Remember to eliminate rows for which the subdeterminant is negative.


#PARAMETERS FOR A LARGE DIRTY INPUT TAX (APPENDIX D.3)
eLK<-1
sX<-1
sU<-1
gL<-0.25
gK<-0.25
tZ<-0.5
tL<-0.3
tK<-1-tL-tZ
sL<-gK*tL/(gK*tL+gL*tK)
sK<-1-sL
d<-gL*sL/tL
Y<-d/(1+d)
L<-sL*(1+gL)/(sL*(1+gL)+sK*(1+gK))
TZHAT<-0.10

##Then run the Table2 function


#PARAMETERS ADAPTED FROM FULLERTON AND TA (2019) (APPENDIX D.4)
eLK<-1
sX<-1
sU<-1
gL<-0.012329
gK<-0.030859
tL<-0.138
tK<-0.415
tZ<-1-tL-tK
sL<-gK*tL/(gK*tL+gL*tK)
sK<-1-sL
d<-gL*sL/tL
Y<-d/(1+d)
L<-sL*(1+gL)/(sL*(1+gL)+sK*(1+gK))
TZHAT<-0.10

##Then run the Table2 function
#Different table structure because since gammaL and gammaK are different, pX and Psources are not equivalent 
xtable(cbind(M2$ekz,M2$elz,M2$pX[,1],M2$pY[,1],M2$merel[,1],M2$w[,1],M2$r[,1],M2$merel2[,1],M2$sd[,1]),caption = "pollu",digits=c(1,1,1,2,2,2,2,2,2,2))
xtable(cbind(M2$ekz,M2$elz,M2$pX[,2],M2$pY[,2],M2$merel[,2],M2$w[,2],M2$r[,2],M2$merel2[,2],M2$sd[,2],M2$pX[,3],M2$pY[,3],M2$merel[,3],M2$w[,3],M2$r[,3],M2$merel2[,3],M2$sd[,3]),caption = "sourceANDuse",digits=c(1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2))
#NOTE: Remember to eliminate rows for which the subdeterminant is negative.


###INCIDENCE ON THE USES SIDE (SECTION 4.2)
##Example where PYHAT-PXHAT<0, and developing some intuition
#This is the second example given in the paper, with positive Allen cross-price elasticities of substitution
tL<-0.5
tK<-0.4
tZ<-1-tL-tK
gL<-4
gK<-1
eLK<-1
eLZ<-1
eKZ<-8
sX<-0.5
sU<-1
sL<-gK*tL/(gK*tL+gL*tK)
sK<-1-sL
d<-gL*sL/tL
Y<-d/(1+d)
b<-Y
a<-sL

eLL<-ell(eLK,eLZ,tL,tK,tZ)
eKK<-ekk(eLK,eKZ,tL,tK,tZ)
eZZ<-ezz(eLZ,eKZ,tL,tK,tZ)

eLL
eKK
eZZ
eLL*eKK-(eLK)^2
eLL*eZZ-(eLZ)^2
eKK*eZZ-(eKZ)^2

SIGN<-sX*(1+gL*sL+gK*sK)+(tL*gK*(1+gL)+tK*gL*(1+gK))*eLK+gL*(1+gK)*(sK-tK)*eLZ+gK*(1+gL)*(sL-tL)*eKZ


DELTA<-sX*(1-b*tZ)*(1+gL*sL+gK*sK)+sU*(gL-gK)*(tL-sL+tZ*(b*sL+(1-b)*a))-eLL*tL*gL*(1+gK)*(b*tK+(1-b)*(1-a))-eKK*tK*gK*(1+gL)*(b*tL+(1-b)*a)+eLK*(tL*gK*(1+gL)*(b*tK+(1-b)*(1-a))+tK*gL*(1+gK)*(b*tL+(1-b)*a))
GENSOURCEINCID<-tZ/DELTA*(sU*(gK-gL)+gL*(1+gK)*eLZ-gK*(1+gL)*eKZ)*TZHAT
GENUSEINCID<-tZ/DELTA*(sX*(1+gL*sL+gK*sK)+(tL*gK*(1+gL)+tK*gL*(1+gK))*eLK+gL*(1+gK)*(sK-tK)*eLZ+gK*(1+gL)*(sL-tL)*eKZ)*TZHAT


#since PHAT=0 (numeraire) we have that WHAT=sK(WHAT-RHAT) -Y(PYHAT-PXHAT)=sK*SOURCE-Y*USE.
WHAT<-sK*GENSOURCEINCID-Y*GENUSEINCID
RHAT<-WHAT-GENSOURCEINCID
LYHATKYHAT<-tL*(eLL-eLK)*WHAT+tK*(eLK-eKK)*RHAT+tZ*(eLZ-eKZ)*TZHAT
LXHATKXHAT<-sX*(RHAT-WHAT)
LYHAT<-(LYHATKYHAT*gK+LXHATKXHAT)/(gK-gL)
KYHAT<-(LXHATKXHAT+LYHATKYHAT*gL)/(gK-gL)
LXHAT<--gL*LYHAT
KXHAT<--gK*KYHAT
XHAT<-sL*LXHAT+sK*KXHAT
PXHAT<-sL*WHAT+sK*RHAT
PYHAT<-tL*WHAT+tK*RHAT+tZ*TZHAT



######################################APPENDIX A.6##########################################

##Examples of pollution-enhancing tax
sU<-0
sX<-0
tL<-1/8
tK<-1/2
tZ<-1-tL-tK
gL<-1
gK<-12/13
eLL<--1
eKK<--1.1005
eZZ<--3
sL<-gK*tL/(gK*tL+gL*tK)
sK<-1-sL
d<-gL*sL/tL
eLK<-(tZ^2*eZZ-tL^2*eLL-tK^2*eKK)/(2*tL*tK)
eLZ<-(tK^2*eKK-tL^2*eLL-tZ^2*eZZ)/(2*tL*tZ)
eKZ<-(tL^2*eLL-tK^2*eKK-tZ^2*eZZ)/(2*tK*tZ)
TZHAT<-1

b<-0
a<-sL
D2<-tZ*tK/sL*(eKK*gK+eZZ*(gK-d)^2/gK-2*eKZ*(gK-d))
DELTA<-sX*(1-b*tZ)*(1+gL*sL+gK*sK)+sU*(gL-gK)*(tL-sL+tZ*(b*sL+(1-b)*a))-eLL*tL*gL*(1+gK)*(b*tK+(1-b)*(1-a))-eKK*tK*gK*(1+gL)*(b*tL+(1-b)*a)+eLK*(tL*gK*(1+gL)*(b*tK+(1-b)*(1-a))+tK*gL*(1+gK)*(b*tL+(1-b)*a))
POLLU<-TZHAT/DELTA*((1+d)*tZ*eZZ*sX+D2*sU-tZ*sX*sU-(gL*tK+gK*tL+gL*gK)*(tL*eLK*eLZ+tK*eLK*eKZ+tZ*eLZ*eKZ))

##Another example with sX and sU different from zero
DELTAFUNC<-function(su,sx,tl,tk,gk,elas,ekk,elk){
  tz<-1-tl-tk
  gl<-1
  ell<--elas
  ekk<-elas*ekk
  elk<-elas*elk
  sl<-gk*tl/(gk*tl+gl*tk)
  sk<-1-sl
  d<-gl*sl/tl
  b<-0
  a<-sl
  delta<-sx*(1-b*tz)*(1+gl*sl+gk*sk)+su*(gl-gk)*(tl-sl+tz*(b*sl+(1-b)*a))-ell*tl*gl*(1+gk)*(b*tk+(1-b)*(1-a))-ekk*tk*gk*(1+gl)*(b*tl+(1-b)*a)+elk*(tl*gk*(1+gl)*(b*tk+(1-b)*(1-a))+tk*gl*(1+gk)*(b*tl+(1-b)*a))
  ezz<-(tl^2*ell+tk^2*ekk+2*tl*tk*elk)/tz^2
  return(cbind(delta,ezz,ell,ekk,elk))
}

DELTAFUNC(0.25,0.25,0.51,0.22,0.128,11,-36,-5.9999)

#Now calculating the change in pZ relative to the CPI and PPI indices in the FH model
#I use the normalization Phat=0 where P=p_X is the FH numeraire, thus pXhat=0
#Therefore the difference pZhat-CPIhat is equal to tauZ-thetaYpYhat
#pYhat can be inferred directly from the uses formula
Y<-d/(1+d)
USE<-tZ/DELTA*(sX*(1+gL*sL+gK*sK)+(tL*gK*(1+gL)+tK*gL*(1+gK))*eLK+gL*(1+gK)*(sK-tK)*eLZ+gK*(1+gL)*(sL-tL)*eKZ)*TZHAT
DIFFCPI<-TZHAT-Y*USE
#Now with the PPI, we know that thetaXL WHAT+thetaXK RHAT=PXHAT=0. We want thetaL WHAT + thetaK RHAT. We have the sources formula which gives WHAT-RHAT
#It is easy to show that thetaL WHAT + thetaK RHAT = (thetaL*thetaXK-thetaK*thetaXL)*(What-Rhat)
#and that thetaL*thetaXK-thetaK*thetaXL = thetaXL*thetaXK(gammaL-gammaK)/(thetaXL(1+gammaL)+thetaXK(1+gammaK))
SOURCE<-tZ/DELTA*(sU*(gK-gL)+gL*(1+gK)*eLZ-gK*(1+gL)*eKZ)*TZHAT
DIFFPPI<-TZHAT-sL*sK*(gL-gK)/(sL*(1+gL)+sK*(1+gK))*SOURCE


##Footnote 14 in F&H 2007
##F&H 2007 assume eLK>0, therefore if ZHAT>0 the numerator must have the wrong sign. Checking by computing D1, D2, D3, D4 and the numerator value
eLK<-2
eLZ<-5
eKZ<--1
sX<-1
sU<-1
gL<-9
gK<-4
tZ<-0.1
tL<-0.18
tK<-1-tL-tZ
sL<-gK*tL/(gK*tL+gL*tK)
sK<-1-sL
d<-gL*sL/tL
Y<-d/(1+d)
L<-sL*(1+gL)/(sL*(1+gL)+sK*(1+gK))

TZHAT<-0.10

eLL<-ell(eLK,eLZ,tL,tK,tZ)
eKK<-ekk(eLK,eKZ,tL,tK,tZ)
eZZ<-ezz(eLZ,eKZ,tL,tK,tZ)

D1<-(1+d)*tZ*eZZ
D2<-tZ*tK/sL*(eKK*gK+eZZ*(gK-d)^2/gK-2*eKZ*(gK-d))
D3<--tZ
D4<--(gL*tK+gK*tL+gL*gK)*(tL*eLK*eLZ+tK*eLK*eKZ+tZ*eLZ*eKZ)
DELTAZ<-D1*sX+D2*sU+D3*sX*sU+D4

b<-0
a<-sL
DELTA<-sX*(1-b*tZ)*(1+gL*sL+gK*sK)+sU*(gL-gK)*(tL-sL+tZ*(b*sL+(1-b)*a))-eLL*tL*gL*(1+gK)*(b*tK+(1-b)*(1-a))-eKK*tK*gK*(1+gL)*(b*tL+(1-b)*a)+eLK*(tL*gK*(1+gL)*(b*tK+(1-b)*(1-a))+tK*gL*(1+gK)*(b*tL+(1-b)*a))

EFFECT<-DELTAZ/DELTA*TZHAT

eLL*eKK-(eLK)^2
eLL*eZZ-(eLZ)^2
eKK*eZZ-(eKZ)^2

